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Abstract 

We extensively develop an algorithm of implementing the Hartree-Fock-Bogolyubov 
calculations, in which the Gaussian expansion method is employed. This algorithm 
is advantageous in describing the energy-dependent exponential and oscillatory 
asymptotics of the quasiparticle wave functions at large r, and in handling vari- 
ous effective interactions including those with finite ranges. We apply the present 
method to the oxygen isotopes with the Gogny interaction, keeping the spherical 
symmetry. In respect to the new magic numbers, effects of the pair correlation on 
the = 16 and 32 nuclei are investigated. 
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1 Introduction 



Experimental data using the radioactive beams have disclosed exotic natures of atomic 
nuclei far from the /3-stability. Some of the data cast questions on the conventional nuclear 
structure theories, which have been founded mainly on top of the data along the /^-stability 
line. Density profile of nuclei near the proton or neutron drip-line is remarkably different 
from that of the stable nuclei, sometimes making halos Furthermore, new magic 
numbers, e.g. N = IQ and 32, have been reported in the region off the /5-stability line 
A property of the effective interaction, which has not been investigated sufficiently, seems 
to play a role in the variation of magic numbers [21 E] • 

The mean-field theories provide us with a good first approximation for bound states 
of nuclei. The Hartree-Fock (HF) theory, through which we can construct the nuclear 
mean-field in a self-consistent manner, is a desirable tool to obtain the single-particle 
(s.p.) orbits of nuclei from microscopic standpoints. In addition, the HF theory is useful 
in investigating basic and global natures of effective interactions. However, the pair corre- 
lation between identical nucleons yields a significant correction to the HF approximation. 
The pair correlation, which gives rise to the superfluidity, can be treated self-consistently 
in the Hartree-Fock-Bogolyubov (HFB) theory For nuclei in the vicinity of the drip 
line, appropriate treatment of the wave-function asymptotics at large r may be important, 
as has been shown in the HF framework [Hj. If the pairing is taken into consideration, the 
s.p. levels that are far below the Fermi energy in the HF scheme couple to the continuum. 
In applying the HFB approximation to nuclei near the drip line, we need a method capable 
of handling the wave-function asymptotics for bound nucleons and those for nucleons in 
the continuum simultaneously. While this point has been argued in association with the 
zero-range Skyrme interaction [Zj, it has not been considered sufficiently for finite-range 
interactions, primarily because no efficient methods to treat the wave-function asymp- 
totics have been known in the case that the non-local nuclear currents are required. 

We recently proposed a new computational method to implement the HF cal- 
culations, by extensively applying the Gaussian expansion method (GEM) [S]. It was 
demonstrated that the energy-dependent exponential asymptotics of wave functions are 
handled efficiently, even with finite-range interactions. In this article we further extend 
the method to the HFB calculations. Both the exponential and the oscillatory asymp- 
totics can be practically represented in a single set of Gaussian bases. This method is 
easily adaptable to calculations with finite-range interactions, and therefore will be useful 
in studying effective interactions to describe nuclear properties including the pair corre- 
lation. The present method is applied to the oxygen isotopes in practice, keeping the 
spherical symmetry. The interaction-dependence of the shell structure near = 16 and 
32 is also investigated, in the light of the pair correlation. 

2 Single-particle bases and HFB equation 

In this paper we assume the mean fields to be spherically symmetric and to preserve the 
parity. Though almost straightforward, extension to the deformed cases will be left as a 
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future study. 

To represent the s.p. wave functions, we introduce bases having the following form: 



Ruij{r) = Kejr^exp{-ur^) . (1) 

Here Y^^\r) expresses the spherical harmonics and Xrr the spin wave function. We drop 
the isospin index without confusion. The index u corresponds to the range parameter 
of the Gaussian basis. We allow the range parameter u to be complex [H], and hereafter 
denote Re(i^) and Im(i^) by u,. and u^, respectively. Note that u,. > 0. The constant f/^ij 
is determined by 

2^+3 21+3 

7ry(2£+l)!! 

so as for {^vijm\^vijnii to be unity. The bases of Eq. (|T)) with different z/'s are not 
orthogonal to one another. The norm matrix for each (^, j) is composed of the elements 
of ^ 

Ar(^i) / I \ /^2yzv^y^2 

Kv' = {^uijmWu'ijm) = I , ^ 1 ■ (3) 

In the GEM we take u's belonging to a geometric progression. The GEM basis-set with 
real u's has been shown to work efficiently in the HF calculations [H] as well as in solving 
few-body problems [TU], including the cases of loosely bound systems. The exponential 
decrease of density at large r is described to a good approximation by a superposition of 
the Gaussians with various ranges. For a complex u, the linear combinations of R^ej and 
Ri,*ej clarify the oscillating structure as 

[R,yej{r) + R^*ij{r)]/2 = J\fui>j / exp(-z/i.r^) cos(i/ir^) , 
[R^tj{r) - R^*tj{r)]/2i = Af^ij exp{-u,r^) sm{iyir^) . (4) 

Indeed, by superposing the complex-range Gaussian bases, oscillating behavior of the wave 
functions in the scattering states as well as in the states with high-nodal wave functions 
can be expressed efficiently [H]- 

To formulate the HFB theory for non-orthogonal bases, we consider the creation (an- 
nihilation) operator cl^j^ {cyijm)-, which is associated with i^utjm and obeys the non- 
canonical commutation relations. 

The HFB equation for non-orthogonal bases is derived generally in Appendix and 
it is reduced to the equation with the i and j conservation in Appendix El Under the 
spherical symmetry, which implies that i and j as well as m are conserved, the generalized 
Bogolyubov transformation is written as 



a 



nljm ^ y 
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where the s.p. state {uijm) stands for the time reversal to (uijm). The solution of 
the HFB equation (PT|) gives the quasiparticle (q.p.) state represented by (nijm), whose 
eigenvalue, enej, is called quasiparticle energy. In Z or = odd nuclei, we should replace 
(U(^J), -V(^J)) by (V(^J)*, U(^^>) in Eq. ^ for a nucleon which has the lowest €nij. With 
the spherical symmetry, this implies addition of 2j+li^un'^^i^n ~^ ^vn^*^v'n) to p"^l) and 

-2kl^^'^^^^'^'n*^^yn^*'^vi) to 4?'^ in Eq. With this modification, the form of the 
HFB equation ()4ip does not change. If we use the linear combination of Eq. instead 
of the complex-range Gaussian bases, all matrices in Eq. (|4ip can be taken to be real. 



3 Asymptotic behavior of quasiparticle wave func- 
tions 

We next consider the asymptotic behavior of the q.p. wave functions at large r. The HFB 
equation (jlT|l is easily represented in terms of the radial coordinate r, by converting the 
label of the basis v to r. For the time being we consider the HFB equation for neutrons. 
At sufficiently large r, the nuclear force becomes negligible, and the following asymptotic 
equations are derived [TT] . 
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'2Mdr^ 



-Aj [rf/l^^-)(r)] ^ 6„,,[rf/l^^-)(r)], 

->] [rV^''\r)] ^ -e„,,[rV;(^^-)(r)]. (7) 



As long as the nucleus is bound, the chemical potential A should be negative. Since we 
take Cnij to be positive, Eq. ^ derives the asymptotic form, 

rV^^^^(r)^ex^(-7i^^'^r) rU^^^^(r) ^ I exp(-r/i^^V) for A + e„^j < 



where rj^^ = ^2M{—X ± e„^j), p^^^ = y^2M(A + e„£j) and 6^^^ is an appropriate real 

number. Note that 77^+'' > rjl^!!} . Corresponding to the asymptotic behavior of Uj^^\r), 
the q.p. energies are discrete only for < e„£j < —A. Thus the deeply bound s.p. states 
in the HF approximation are embedded in the continuum in the HFB theory. It is now 
obvious that, in the HFB calculations, we have to treat both the exponential and the 
oscillatory asymptotics which depend on the q.p. energy. In Ref. [7], these asymptotics 
are handled by adopting the transformed harmonic oscillator (THO) bases. However, it 
is not easy to apply the THO bases if the non-local currents are required, as in the case 
of the finite-range interactions. As will be demonstrated in Sec. the GEM provides a 
practical method in handling the energy- dependent asymptotics even in the presence of 
the non-local currents. 

The asymptotic behavior of the neutron density and pair current jT^j is derived from 
Eq. (jS)). The density matrix and the pairing tensor in Eq. (|39p are converted to the 
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r-representation as 



n 

^m^ry) = Y.V^'^^*ir)Ui'^\r'), (9) 

n 

for = even nuclei. The density and the pair current are then defined by 

p(r) = ^(2j + l)p(^^)(r,r), 

K{r) = ^(2j + l)fi:(^J)(r,r). (10) 
ij 

Denoting the smallest q.p. energy by e^i^, we immediately obtain the asymptotic form of 
p(r) as 

rV(r) ^ exp{-27]';^'''r) , (11) 



where 77™™ = y 2M(— A + emin)- In order for an = odd nucleus to be bound, the 
condition A + emin < must be fulfilled, because Vj;^^^\r) is replaced by Uj{^\r) for the 
last neutron. The contribution of the last neutron yields the asymptotic form 

rV(r) ^ exp(-2r/!;^'°r) , (12) 

which damps more slowly than the form in Eq. (fTT|) . 

As pointed out in Ref. j^, a complication may arise in the close vicinity of the drip 
line. If —A is almost vanishing and no q.p. level satisfies X + enej < 0, it is difficult to draw 
definite conclusion on the asymptotic behavior of p(r). For A^ = even nuclei, arguments 
based on the q.p. energy suggest that the component of X + enij +0 becomes dominant 

However, 



at the r — > 00 limit, giving the asymptotics of r"^ p{r) ^ exp — 2y4M(— A) r 
the spectroscopic amplitude of this component may be negligibly small in the physically 
interesting region. The actual situation is in the balance between the damping factor 
and the spectroscopic amplitude of the q.p. levels near —A. Moreover, because of this 
subtlety, in practical calculations the asymptotic behavior of p(r) could be sensitive to 
the details of the method, even if individual q.p. wave function has the asymptotics 
consistent with its q.p. energy. Even with the high adaptability of the GEM to the wave- 
function asymptotics at large r, it is beyond scope of the present study to obtain reliable 
asymptotics of p(r) for nuclei with A + emin > 0. Similar arguments apply to A^ = odd 
nuclei. 

The asymptotics for K{r) could also be complicated. We separately consider contribu- 
tion of the discrete q.p. states and that of the q.p. continuum. The discrete q.p. states, 
which satisfy A + e„£j < 0, contribute to n{r) asymptotically with the damping factor 



exp 



{f]n+ + Vn-)^ ■ Because rj^l' + r)ni' is a decreasing function of enij for a fixed A, 
the q.p. level just below —A dominates over the other discrete states at large r. The 



damping factor becomes exp — 2i/4M(— A)r at the A + enij —0 limit. In the states 
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having A + e„£j > 0, we have a damping factor exp(— ?7^+V), which is muhiphed by the 
oscillating factor cos(p^^-'V + 6'^^-''''), at large r. Hence the q.p. level closest to —A has dom- 
inant contribution to K(r) at the r — > oo limit. The component of A + e„£j — > +0 gives 



the most slowly damping factor exp — y4M(— A)r . However, similarly to p(r) in the 
A ~ case, spectroscopic amplitudes of the q.p. levels become relevant. This complicates 
description of the asymptotic behavior of K{r) near the neutron drip line. 

For proton wave functions, the Coulomb interaction influences the wave functions at 
large r, although it does not affect the criterion whether or not individual q.p. levels are 
discrete. The above arguments for the neutron functions can be applied if the asymptotic 
forms are properly modified. Nevertheless, the complication near the drip line is expected 
to be less serious than for neutron wave functions, owing to the presence of the Coulomb 
barrier. 



4 Effective interaction 

We shall consider the effective Hamiltonian comprised of the kinetic energy and the ef- 
fective two-body interaction, 

H = K + V- K = Y.^, V^ = E%- (13) 

i ^^^^ i<j 

Here i and j are the indices of each nucleon. The s.p. matrix element of the kinetic term 
is calculated as 

As in Refs. jH E] , the effective interaction Vij is assumed to have the following form. 



{V^uijm 1 7777 1 ^-'ijrn) = ■ -r—; KJ' ■ (14) 



^12 — t^l2 ^ ^12 ^ ^12 ^ ^12 ) 
^12 



„,(LS) 
'^12 



= E(f "^^SE + C^PtE + tf^^PsO + tr'PT0)fP{n,) , 

= E(^f "^^TE + ^°)PT0)/r(--) . (S, + S,) , 
^(TN) ^ 5:(t(™E)p^^ + ,(TNO)p^^)^(TN)(^^^)^2^^^^^ 

= t,il + xsP.Mr,)r5{vu), (15) 

where = ri - r2, ri2 = |ri2|, P12 = (pi - P2)/2, L12 = ri2 x pi2, Sj is the nucleon 
spin operator, and 5*12 is the tensor operator Su = 4 [3(si ■ fi2)(s2 ■ 1-12) — Si ■ S2]. The 
projectors Pse, Pte, Pso and Pto are expressed as 

^ _ 1 - 1 + „ _ 1 + P, 1 - P, 

PSE-^ 

I- Pa I- Pr l+Pa l + Pr 

Pso = , Pto = , (16) 
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in terms of the spin and isospin exchange operators P^, Pt- In Eq- (US)), //^ represents 
an appropriate function, /x stands for the parameter attached to the function, and the 
coefficient. By taking various functions for Eq. ()15|1 covers a wide variety of effective 
interactions. The Gogny interaction jTHj is obtained by setting f^\ri2) = exp(— /xr^g); 
f^^^\r 12) = V^5(ri2) and fj™'' = 0. In Ref. jlj we have developed M3Y-type interactions, 
in which we take /^^'^^'™'*(ri2) = exp(— yuri2)//iri2. 

The matrix elements of w^2 '^^'™'' can be calculated by utilizing the Fourier transfor- 
mation of ffj,{r 12) [HI- Most formulae in Ref. [H] are straightforwardly applicable, although 
the complex-range Gaussian bases lead to a slight modification in the radial integral 
(Eqs. (24-26) in Ref. 0): 



00 

2 







r drjx{kr)R^* 



2 2 - 




'(2£+ l)!!(2f + 1)!! {u* + u')'^ 
where LI^\x) is the associated Laguerre polynomial and 

Q 



■ Ql., L^'-.ll (Ql,,) exp{-Ql.,,) , (17) 



2(u* + u')^ 



We need only the i + f + \ = even cases for the integral of Eq. (fTTjl . In coding a computer 
program, we prepare a subprogram for the fc- integration (see Eq. (28) of Ref. [H]), which 
is the only part dependent on the form of /^(t'i2)- This allows us to handle various 
interaction forms only by substituting the subprogram. Although analytic formulae for 
the A;- integration were given in Ref. jH] , their application sometimes causes serious round- 
off errors. Numerical integration over k is recommended to reduce the errors. 

Mean-field approaches for nuclei are sometimes discussed in connection to the density- 
functional theory (DFT). Although the DFT is a rigorous theory, we have to know a good 
density-functional form of energy for practical applications. The most popular density- 
functional form in nuclear structure physics is based on the Skyrme interaction and its 
extensions, whose parameters are adjusted to the known data. Some of the Skyrme density 
functionals involve currents that cannot be derived from two-body interactions. Still the 
form of the Skyrme density functionals is quite restricted, since all currents are assumed 
to be local. In order to know appropriate form of density functional for nuclei covering 
large area of the nuclear chart, it will be an important step to explore finite-range effective 
two-body interactions. It is noted that, whereas we consider two-body interactions having 
the form of Eq. (jl5|) in this paper, the GEM algorithm is applicable also to the DFT, 
which may not necessarily be connected to two-body interactions. 

Once storing the interaction matrix elements and having a trial set of f/^^^^ and Vj!^^\ 
we obtain the HFB equation pij) with /i^f^'' and A^f^^ defined in Eqs. and (PUj). and 
the equation is solved iteratively until convergence. 
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5 Numerical examples 



We now apply the present method of the HFB calculations to several nuclei, assuming 
the spherical symmetry and the parity conservation. The center-of-mass energy is fully 
removed before variation, by subtracting both the one-body and two-body terms from the 
Hamiltonian of Eq. ()13|) . 

5.1 Application to oxygen isotopes with Gogny interaction 

We first show results for oxygen isotopes, using the standard DIS parameter-set [Hj of the 
Gogny interaction. While the real-range GEM basis-sets are efficient to describe nodeless 
but broadly distributing wave functions, the complex-range GEM sets seem suitable for 
reproducing the oscillating behavior required in reaction studies jHj. In the HFB calcu- 
lations both types of wave functions come into the problem. We try the following three 
basis-sets, a set comprised of real-range bases, a set of complex-range bases and a set of 
their mixture; 

Set A: = ly^b"^'^ , 6 = 1.20 
SetB: = 1.^6-2", 6 = 1.15 

SetC: = z/^6-2", 6=1.25 

where z/^ = Mui /2h with hw = 41.2 x 24"^/^ MeV. The parameter K stands for number 
of bases for each (£, j), and is fixed to be 12 for all the three sets. We truncate the 
space, taking the £ < 4 s.p. bases. If the £ = 5 bases are added, we have energy gain of 
< 0.06 MeV. 

In practical calculations, the common ratio 6 of the GEM bases is restricted so that the 
norm matrix should not have a singularity. At a certain value of 6 which is close to unity, 
convergence in the iterative process becomes extremely slow in numerical calculations with 
a given precision. We here call this phenomenon 'quasi-instability'. If 6 further approaches 
unity, true instability comes out, in which round-off errors lead to unphysical results, 
because the norm matrix becomes nearly singular. Presently the parameter 6 is chosen 
for each set so as to optimize the energy of ^^O, with avoiding the quasi-instability in 
calculations using the double precision. In the present calculation ^^O is bound although it 
lies beside the neutron drip line, whereas the actual ^^O seems unbound ^3]. By contrast, 
2^0 is unbound in the respect that the neutron chemical potential becomes positive, which 
leads to the nuclear wave function distributing up to the infinite distance. For Set B, the 
value of Vi/v-r has also been tuned. In the GEM algorithm of the mean-field calculations, 
the parameters for the basis-set do not need to depend strongly on mass number. Using 
the same basis-set, we can obtain results for all nuclides under consideration in a single 
run, without recalculating the matrix elements of the effective interaction. 

If the DIS force is applied to the pure neutron matter, the energy per nucleon diverges 
with the negative sign at the high density limit. Due to this nature, the true energy 



z^i = , (a 
^ = ±-■0.6, (a 



-2,-l,0,--- ,J^-3); 
0,1,--- ,K/2-l); 

0,---,K/2-l) 
O,--- ,ir/4-l) ' 
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minimum in a finite nucleus lies at the configuration in which all neutrons gather in the 
vicinity of the origin {i.e. the center-of-mass) without overlap of the proton distribution. 
To avoid the tunneling to this unphysical configuration, the range of the GEM bases 
(z/jT^/^) should not be too small. All the present sets are chosen to meet this criterion. 

We tabulate total energies of i^-^gq obtained from the HFB calculations using the 
three basis-sets, in Table ^ Energies in the HP approximation are also presented. The 
result for ^^O of the HFB calculation with Set A is not shown, because the quasi-instability 
takes place. The proton Os- and Op-shells are fully occupied, forming the Z = 8 closed core, 
in any results of Tabled In i5-i7,2i,23-25q^ neutrons lose the superfluidity, irrespectively 
of the basis-sets. Because the variational principle is satisfied in the HP and HPB theories, 
the lower total energy indicates the better result for each nucleus. The energies obtained 
from Sets A and B are comparable, with difference less than lOkeV in most oxygen 
isotopes. To be more precise, while it depends on the nuclide which of Sets A and B gives 
lower energy in the HP calculation. Set B yields steadily lower energy than Set A in the 
HPB calculation, with the exception of ^^O. This is obviously because the energy gain due 
to the pairing is larger in Set B than in Set A. Giving energies lower than Sets A and B by 
several tens keV, Set C is advantageous over the other sets for any oxygen isotopes in both 
the HP and the HPB calculations. Because the energy difference among the basis-sets is 
contributed by individual q.p. or s.p. orbits, difference in total energy will roughly be 
proportional to the mass number, growing significantly for heavier nuclei fB] • 

We next look into the wave functions of the q.p. states, using Set C. Figure [T] shows 
the radial wave functions of the neutron q.p. levels in ^^O which correspond to the Osi/2, 
lsi/2 and 2si/2 orbits, in terms of \rV^^^\r)\'^ and \r Uj{'^\r)\'^ . In this calculation ^^O 
is so loosely bound with A ^ — 0.4MeV that all q.p. levels would lie in the continuum, 
satisfying A -|- enej > 0. In this regard the q.p. levels shown in Pig. Q are considered 
to be resonance- like levels. For all the three q.p. states, whose q.p. energies distribute 
from 1.4 MeV to 40MeV, the energy-dependent exponential asymptotics in rV^^^\r) are 
represented appropriately by the superposition of the GEM bases. We find oscillatory 
behavior in r U^^\r). In order to view how appropriately the oscillatory asymptotics are 
described, we define the following quantity from the Fourier transform of rUj^^\r), 

' {=T^eA-k)) ■ (19) 

The scale of Tn£j{k) is connected to the occupation probabihty of the q.p. level, via 

/ TnN{k)dk= \Ui'^\r)\^r^dr = l- \V^'^\r)\^ r^dr . (20) 
JO Jo Jo 

If Uj{^\r) fulfills the asymptotics of Eq. (jHl) with A + e„^j > 0, Tnij{k) should have a 
peak at k = pl^^^ . In Fig. |21 we show Tnej{k) for several low-lying q.p. levels of ^^O. 
Viewing the peaks of Tnej{k) at the positions compatible with the q.p. energies, we 
confirm that the energy-dependent oscillatory asymptotics are properly expressed through 
the superposition of the GEM bases. 



^nljik) = — 
TX 
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In Table |21 we show the rms matter radii in HF and HFB obtained through the three 
basis-sets. Contribution of the center-of-mass motion is subtracted as in Ref. There 
is no significant difference among the bases up to ^^O, and this seems consistent with 
the results in Tabled However, we find sizable basis-dependence in 25.26q^ although the 
energies do not differ significantly. In ^^O, where neutrons are not in the superfluid phase, 
the s.p. energy of n0d3/2 in the HF scheme is vanishing — O.OSMeV). Thereby the 
calculated radius of ^^O is extremely sensitive to precise value of the s.p. energy. The ^^O 
case gives an illustration of the subtlety for the drip-line nuclei that have no discrete q.p. 
levels, as has been discussed in Sec. El In these cases in which the nuclei are very close 
to the neutron drip line, some physical quantities {e.g. nuclear radius) may depend on 
details of the calculation. Comparison of the radii between the HF and HFB calculations 
will be of a certain interest. The matter radii in the HFB scheme are similar to those in 
HF scheme around the /^-stability line. However, in ^^O the HFB calculations give slightly 
larger radius than the HF calculations, primarily due to the partial occupation of nlsi/2- 
Contrastingly, in ^^O the HFB radius is substantially smaller than the HF one in Sets B 
and C. This corresponds to the pairing anti-halo effect ^2), although attention should be 
paid to the subtlety discussed in Sec. El 

5.2 Dependence of shell structure on HF interaction 

It has been argued that the spin-isospin nature of the A^A^ interaction could play a role 
in producing the new magic numbers [3]. We have developed an M3Y-type effective in- 
teraction called M3Y-P2 jl], which is applicable to the self-consistent HF calculations. 
Keeping the one-pion-exchange contribution to the central part of the NN interaction, 
the M3Y-P2 interaction possesses reasonable spin-isospin nature, as was checked via the 
Landau-Migdal parameter in the nuclear matter. Within the self-consistent HF calcu- 
lations, shell gaps at = 16 |3] and at = 32 ^Tj behave quite differently between 
M3Y-P2 and frequently used interactions such as DIS. With the DIS interaction the shell 
gap between nlsi/2 and nOds/2 in the A^ = 16 isotones is almost constant, as Z changes 
from 8 to 16. On the contrary, the shell gap becomes narrower as Z increases, if we adopt 
the M3Y-P2 interaction. Similar behavior has been found in the A^ = 32 isotones. It 
should be pointed out that major difference lies around the /3-stable region, rather than 
in the region close to the drip-line. We confirmed that the Z-dependence in the M3Y-P2 
result is primarily ascribed to the contribution of the one-pion-exchange part. 

The pair correlation can break the closed core, even in the spherical nuclei. This core 
breaking is realized in the HFB picture, unless the relevant shell gap is large enough. To 
investigate effects of the pair correlation in the region concerning the new magic numbers, 
we carry out the HFB calculations using the GEM algorithm for the A^ = 16 and 32 
nuclei. While the DIS interaction has moderate pairing properties, the pairing character 
of the M3Y-P2 interaction is not realistic, since it has not been tuned. Here we always 
use the DIS interaction for the pairing interaction, i.e. to obtain the pair potential in the 
HFB equation ()41|). and examine dependence of the magic nature of A^ = 16 and 32 on 
the HF interactions, for which we take DIS and M3Y-P2. The ^ < 4 < 5) truncation is 



9 



made for the = 16 (A^ = 32) nuclei. For the basis-set, Set C in the preceding subsection 
is employed, with K = 12. 

In TableEl the neutron pair energies in the = 16 and 32 nuclei are compared between 
the two HF interactions. The pair energy is defined as contribution of the pairing tensor 
to the energy of the system. In the iV = 16 isotones, the DIS interaction does not give 
rise to the superfiuidity for neutrons in 8 < Z < 14, as long as the spherical symmetry 
is assumed. In this regard, the DIS interaction provides a prediction that the = 16 
magic number holds in the whole region of 8 < Z < 14. On the contrary, in the case that 
we use the M3Y-P2 interaction for the HF potential, the neutron pair energy increases as 
Z goes from 8 to 14, i.e. from the drip-line region to the /9-stable region. This behavior 
of the neutron pair energies, which is a refiection of shell structure in the HF scheme, 
is typically viewed in the difference between ^^O and '^'^Si. Similar trend is seen in the 
A^ = 32 isotones. While in ^^Ca the pair energy is not large and is close to each other 
between the two HF interactions, we find sizable interaction-dependence in ^'^Ni, due to 
the difference in the shell structure in the HF regime. 

Effects of the pair correlation on excitation spectra are basically viewed in the q.p. 
energies. In the upper panel of Fig. 01 the neutron q.p. energies corresponding to the 
lsi/2 and 0(1^/2 orbits are compared between the DIS and the M3Y-P2 interactions for 
the Z = even, A^ = 16 isotones. For reference, we also plot the s.p. energies measured 
from the Fermi energy in the HF approximation without the sign, as counterparts to 
the q.p. energies. We here define the Fermi energy as the arithmetic average of the HF 
energies of nlsi/2 and n0d3/2. Hence, half the shell gap between nlsi/2 and n0d3/2 is 
presented as the HF results. In the case that the neutron pair energy vanishes in the 
HFB calculations, we take the q.p. energies to be equal to their HF counterparts. We 
also show the values obtained from a HF picture within the srf-shell, employing the 'USD' 
shell-model Hamiltonian ^Hl- As pointed out in Ref. jH], the HF results indicate that 
the USD interaction yields significant Z-dependence of the shell gap, having qualitative 
similarity to the M3Y-P2 interaction. It should be noticed that the HF treatment of the 
shell-model Hamiltonian tends to give radical variation of the s.p. energies from nuclide 
to nuclide, in comparison with the self-consistent mean-field approaches, because effects 
of the rearrangement of the core are incorporated into the interaction within the valence 
shell in an effective manner. 

The HF and the HFB energies coincide up to Z = 14 in the calculations using the DIS 
interaction, because the nuclei are not in the superfluid phase. In the HFB results, the 
M3Y-P2 interaction provides almost constant q.p. energies of nlsi/2 and n0d3/2, which 
are close to those given by the DIS interaction. This is because the narrowing of the shell 
gap is compensated by the pairing gap to a considerable extent. It is suggested from the 
q.p. energies that the two interactions yield similar excitation spectra. In other words, 
the M3Y-P2 interaction does not seem to destroy the energy spectra predicted in the 
widely used interactions like DIS, despite the difference in physical contents. However, 
additional correlations could affect the energy levels considerably. It will be interesting 
to investigate influence of the correlations beyond HFB in these nuclei. 

The occupation numbers of the relevant q.p. levels are depicted in the lower panel 
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of Fig. El Due to the narrowing shell gap, we observe depletion of the nlsi/2 occupation 
in the HFB results with M3Y-P2, as Z grows from 8 to 14. This indicates, consistently 
with the behavior of the pair energies, that the magic nature of = 16 becomes weaker 
as Z increases. Thus the spectroscopic factors of the neutron orbits could be useful in 
discriminating the two pictures; a picture based on the almost constant iV = 16 shell 
gap, and that on the decreasing shell gap for increasing Z which may be compensated by 
the pairing gap. As shown in Fig. El the occupation numbers obtained with M3Y-P2 are 
compared to those in the full shell-model calculation using the USD Hamiltonian, except 
at Z = 12 where influence of the quadrupole deformation seems important. 

Figure m shows the q.p. energies and the occupation numbers of the levels correspond- 
ing to the nlp3/2, ^0/5/2 and nlpi/2 orbits, for the Z = even, = 32 isotones. As the 
HF counterparts to the q.p. energies, we plot the s.p. energies measured from the Fermi 
energy, without the sign. The Fermi energy is taken to be the arithmetic average of the 
energies of the two levels nearest the = 32 gap {i.e. 1^3/2 and the lower level out of 
0/5/2 and 1^1/2)- As a result, the two levels have equal energy, which represent half the 
shell gap, in regard to the HF counterparts to the q.p. energies. These values are shown 
by the thin black lines in the upper panel of Fig. |3 As argued in Ref. f7j, the DIS in- 
teraction gives almost constant shell gap in 20 < Z < 28, while M3Y-P2 yields erosion of 
the gap as Z goes from 20 to 28. Furthermore, the s.p. energy difference between nO/5/2 
and nlp3/2 changes more rapidly in the M3Y-P2 results than in the DIS ones. This gives 
rise to the inversion of nO/5/2 and nlpi/2 near Z = 24 in the M3Y-P2 case, while such 
inversion does not occur in the DIS case. This shell structure in the HF scheme affects the 
HFB occupation numbers. Although the = 32 core does not close in the whole region 
of 20 < Z < 28 in the strict sense of the HFB theory, at Z = 20 the occupation number 
of each orbit is not significantly different from the value assuming the A^ = 32 closure, 
irrespectively of the HF potential. If we use the DIS interaction for the HF potential, 
this nature is almost preserved as Z increases. In contrast, when we apply the M3Y-P2 
interaction, leakage out of 71.1^3/2 appreciably grows for increasing Z, mainly due to the 
pair excitation to riO/5/2. Despite this difference in the physical picture, the q.p. ener- 
gies are much less interaction-dependent, because the decreasing shell gap in M3Y-P2 is 
compensated by the increasing pairing gap to a certain extent. Therefore the excitation 
spectra will not differ much between the two HF interactions, unless disturbance due to 
the additional correlations is substantially strong. Yet the inversion of nO/5/2 and nlpi/2 
occurs in the HFB results with M3Y-P2, which is not seen in the results with DIS. 

6 Summary 

Extending the previous work on the Hartree-Fock (HF) calculations based on the Gaus- 
sian expansion method (GEM), we develop a new method of implementing the Hartree- 
Fock-Bogolyubov (HFB) calculations. This method has two notable advantages: (i) It is 
efficient in describing the exponential and the oscillatory asymptotics of wave functions 
at large r, which depend on quasiparticle energy and could play a significant role in nuclei 
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close to the drip lines, (ii) We can handle various effective interactions, including those 
inducing non-local nuclear currents. 

The present method has been tested numerically in the oxygen isotopes with the Gogny 
DIS force, by assuming the spherical symmetry and the parity conservation. Three sets 
of GEM bases are tried and their results are compared; a set comprised of the real-range 
Gaussians, a set comprised of the complex-range Gaussians, and a set of their mixture. 
The last one gives the best results in the context of the variational principle. We have 
confirmed that both the exponential and the oscillatory asymptotics of the quasiparticle 
wave functions are properly described. 

Relating to the new magic numbers in nuclei far from the /3-stability, we have investi- 
gated interplay of the pair correlation and the shell gaps at = 16 and 32. The Gogny 
DIS and the M3Y-P2 interactions are used for the HF interaction, while for the pair 
potential we always use the DIS interaction. Reflecting the interaction-dependence of the 
shell structure in the HF scheme, the M3Y-P2 case provides stronger Z-dependence of the 
superfluidity at = 16 and 32 than the DIS case. While the q.p. energies do not differ 
significantly between the DIS and the M3Y-P2 interactions, suggesting that excitation 
spectra could resemble each other, the occupation numbers of neutron orbits reflect the 
interaction-dependence of the shell structure. The occupation numbers obtained from 
the M3Y-P2 interaction are comparable with those in the standard shell model, for the 
TV = 16 nuclei. 

The author is grateful to Y. R. Shimizu for discussions. This work is financially supported 
in part as Grant-in-Aid for Scientific Research (B), No. 15340070, by the Ministry of 
Education, Culture, Sports, Science and Technology, Japan. Numerical calculations are 
performed on HITAC SRllOOO at Institute of Media and Information Technology, Chiba 
University, and at Information Technology Center, University of Tokyo. 

Appendices 

A HFB equation for non-orthogonal bases 

In this appendix we present the HFB theory in the case that non-orthogonal s.p. bases are 
adopted. We denote the creation (annihilation) operator associated with the s.p. basis k 
by c|. (cfc). The non-orthogonality leads to the non-canonical commutation relations. 



where N = {Nkk') is the norm matrix satisfying N^i^, = Nyk and being positive-definite. 
The generalized Bogolyubov transformation is defined by 



{cfc, 4,} = Nkk' , {cfc, Cfc'} = {4, 4'} = , 



(21) 




(22) 



k 
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In the matrix representation, Eq. ()22|) can be expressed as 



(at «) = (ct c)W; \N=(^y ^Ij- (23) 

Obtained as a solution of the HFB equation, a] and ai obey the canonical commutation 
relations, 

{ai, a],} = 6ii' , {ai, ai'} = {a}, a},} = . (24) 



This indicates 



I.e., 



WN'W = 1, WW=N'-^ N'=(q °, J , (25) 



UNU + VN*V = 1, V^NU + U^N*V = 0, 

UU + V*V^ = N~S UV + V*U^ = 0. (26) 

We here use the expression A = (A*)""", instead of A^, to avoid confusion with the hermitian 
conjugation of operators. The inversion of the Bogolyubov transformation of Eq. (j^Hj) 
yields 

( ct c ) = ( at a )WN'. (27) 
We define the density matrix and the pairing tensor by 

(NpN)fcfc, = ($|ct,cfc|$) = (NV*VTN)fefc,, 

= ($|cfc.Cfc|$) = ($144,1$)* = (NV*UTN*)fcfc,, (28) 

where |$) is the HFB vacuum. This HFB vacuum |$) satisfies anijml^) = and ($|$) = 
1. Equation obviously indicates 

Pkk' = (V* V^)kk' , f^kk' = (V* U^)fcA:' . (29) 

In the HFB approximation the expectation value of the number operator is adjusted to 
the actual particle number n. This constraint is expressed by Tr(N p) = n. The HF 
Hamiltonian and the pair potential are defined by 

OPkk' 

^kk' = j^mm) = {mck,H],ck'm, (3o) 

from which the following HFB Hamiltonian matrix is obtained. 
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with the chemical potential A. The HFB equation for the non-orthogonal bases has the 
following form of the generalized eigenvalue equation, 

7^W=N'W(^J ^^y, e = diag(e,). (32) 

It is commented that, because N is a Hermite positive-definite matrix, it can be decom- 
posed as N = L L by using a lower-triangular matrix L. Multiplication of c by gives 
an orthogonalization of the s.p. bases, and the canonical commutation relation follows; 

{(L-ic)fc,((L-i)*ct),,} = 4fc'. (33) 

All the equations return to those in the usual HFB theory for the orthogonal bases, by 
multiplying by appropriately. 

For the Hamiltonian consisting of the one-body term K and the two-body interaction 
V, we obtain 

= {hlKlh) pk,k, + \ {hh\V\kik2) {2pk,k,Pk2k, + t^k,kXk;k,) , (34) 

kik2 k\k2kzk/i 

where \k) = c^|0) and \kk') = c^c^/ 10), with |0) representing the particle vacuum satisfying 
Cfc|0) = 0. Equation derives 

hkk' = {k'\k\k) + 'Y{k'k2\V\kki) pkik2 , 

felfe2 

^fefc' = ^'Y{kk'\V\kik2) Kkik2 ■ (35) 

kik2 

For interactions having the form of Eq. rearrangement terms for v^^\ which emerge 
via variation of p°, should be taken into account. 



B HFB equation with spherical symmetry 

By adopting the spherically symmetric s.p. bases k = {uijm) and k' = {u'i'j'm'), the 
norm matrix N is simplified as in Eq. Q, 

Nkk' = Su'Sjj'^mm' ■ Njfji^ . (36) 

If we assume the spherical symmetry on the HFB fields, the q.p. state is taken to be 
i = {nijm), and the U- and ^-coefficients in Eq. fl22|) have the form 

Uk'i = See'Sjj>6mm' ■ Ul^'n 1 ^k'i = ^U'^jj'^m-m'{ — y^^ ■ Vyin ■ (37) 

Equation is then reduced to 

UM n(^^') u^^^') + \mH^^^> v(^^') = 1 , v(^^')^ n(^^') u^^^') - u^^^')^ n(^^> v(^^') = , 
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Analogously to Eq. ()37|). the density matrix, the pairing tensor, the HF Hamiltonian and 
the pair potential are taken as 

^kk' = Su'S,,'6m-m' i-y-"^ ■ ^i'S ; = (V('^> U(^^)T),,. , (39) 

and 

hkk' = Su'Sjj'^mm' ■ h^ul' ' ^kk' = ^U' ^jj'^m-m' { — Y'"^ " ^l^u' ■ (40) 

Note that nlf^} = k^^/J and A^f^'' = ^u''',!- With these replacements the HFB equation 
under the spherical symmetry is written as 

n^^^) W(^^) = N'(^^) W(^^) [ ^ _ ) ; e(^^) = diag(6„,,) , (41) 



where 
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Table 1: HF and HFB energies (MeV) of ^^"^'^O, calculated with the PIS interaction. 



nuclide 




HF 






HFB 






Set A 


Set B 


Set C 


Set A 


Set B 


Set C 


14Q 


- 99.522 


- 99.533 


- 99.559 




- 99.586 


- 99.613 


150 


-114.467 


-114.479 


-114.509 


-114.467 


-114.479 


-114.509 


160 


-129.477 


-129.483 


-129.515 


-129.477 


-129.483 


-129.515 




-134.518 


-134.521 


-134.548 


-134.518 


-134.521 


-134.548 


180 


-139.711 


-139.710 


-139.734 


-142.283 


-142.286 


-142.309 


190 


-145.079 


-145.076 


-145.096 


-146.391 


-146.391 


-146.411 


20O 


-150.644 


-150.639 


-150.657 


-153.329 


-153.331 


-153.350 


210 


-156.425 


-156.419 


-156.434 


-156.425 


-156.419 


-156.434 


220 


-162.438 


-162.432 


-162.444 


-162.540 


-162.546 


-162.563 


230 


-165.278 


-165.298 


-165.321 


-165.278 


-165.298 


-165.321 


240 


-168.540 


-168.568 


-168.595 


-168.540 


-168.568 


-168.595 


250 


-168.082 


-168.084 


-168.134 


-168.082 


-168.084 


-168.134 


260 


-167.942 


-167.935 


-167.995 


-169.248 


-169.266 


-169.304 



Table 2: Rms radii (fm) of ^^O in the HF and HFB calculations with the DIS interac- 
tion. 



nuclide HF HFB 

Set A Set B Set C Set A Set B Set C 



140 


2, 


.532 


2. 


,530 


2, 


,530 






2. 


,535 


2, 


,534 


150 


2, 


,570 


2. 


,568 


2, 


,568 


2, 


,570 


2, 


,568 


2, 


,568 


160 


2, 


.609 


2. 


,607 


2. 


,606 


2. 


,609 


2. 


,607 


2. 


,606 




2, 


.663 


2. 


,662 


2, 


,661 


2. 


,663 


2. 


,662 


2. 


,661 


180 


2, 


.712 


2. 


,711 


2, 


,711 


2, 


.707 


2, 


.706 


2, 


,705 


190 


2, 


.756 


2. 


,755 


2, 


,755 


2, 


.755 


2, 


.754 


2, 


,753 


20O 


2, 


.796 


2. 


,795 


2, 


,795 


2, 


.794 


2, 


.793 


2, 


,793 


210 


2, 


.832 


2. 


,831 


2, 


,831 


2, 


.832 


2, 


.831 


2, 


,831 


220 


2, 


.864 


2. 


,863 


2, 


,863 


2, 


.875 


2, 


.874 


2, 


,874 


230 


2, 


.950 


2. 


,947 


2, 


,948 


2, 


.950 


2, 


.947 


2, 


,948 


240 


3 


.016 


3. 


,013 


3. 


,014 


3 


.016 


3 


.013 


3. 


,014 


250 


3 


.133 


3. 


,102 


3. 


,119 


3 


.133 


3 


.102 


3. 


,119 


260 


3 


.205 


3. 


,177 


3. 


,197 


3 


.203 


3 


.159 


3. 


,164 



Table 3: Interaction-dependence of neutron pair energies (MeV) in = 16 and 32 nuclei. 



nuclide 


DIS 


M3Y-P2 


240 


0.000 


0.000 


30gi 


0.000 


-3.911 


52Ca 


-2.122 


-2.224 


60Ni 


-4.825 


-7.055 




2 4 6 8 10 12 14 

r [fm] 

Figure 1: Wave functions of neutron Si/2 Q-P- levels in ^^O calculated with the DIS 
interaction. The colors distinguish the q.p. levels: The red, blue and green lines show the 
wave functions of the levels corresponding to the Osi/2, lsi/2 and 2si/2 orbits, respectively. 
The full (dashed) curves represent \r Vji^^\r)\'^ Uj{^\r)\'^) for each level. 




Figure 2: Tnijik) of neutron q.p. levels in ^^O calculated with the DIS interaction. The 
values of p^^-'^ for individual q.p. levels are presented by the ticks on the horizontal axes. 




Figure 3: Neutron quasiparticle energies (upper panel) and neutron occupation numbers 
(lower panel) in the Z = even, = 16 isotones. In both panels the red lines present the 
values of 0c?3/2 in the HFB calculations, while the blue lines those of lsi/2 if visible. The full 
(dashed) lines with each color are obtained with the M3Y-P2 (the DIS) HP interaction. 
The thin black lines indicate the HF results (see text for details). The black dotted line in 
the upper panel is obtained from the HF treatment of the USD Hamiltonian. The dotted 
lines with each color in the lower panel represent the occupation numbers obtained from 
the full shell-model calculation with the USD Hamiltonian. 



3. On 




20 22 24 26 

Z 



Figure 4: Neutron quasiparticle energies (upper panel) and neutron occupation numbers 
(lower panel) in the Z = even, N = 32 isotones. In both panels the red, blue and 
green lines present the values corresponding to the O/5/2, lpi/2 and 1^3/2 levels. The full 
(dashed) lines with each color are obtained with the M3Y-P2 (the DIS) HF interaction. 
The thick lines indicate the HFB results, while the thin lines the HF results. In the upper 
panel, the black color is used for the energy of the two levels adjacent to the Fermi energy 
in the HF results. 



